#####对HMM富集结果meta分析后可视化

library(dplyr)
setwd("E:/5hmc_file/脑组织和GM12878的chromHMM数据")
anno=read.table("group意义.txt",head=T,sep="\t")
anno=data.frame(HMM.group=anno$STATE.NO.,description=anno$DESCRIPTION,num=anno$NUM)

file=read.table("result_HMM_enrichment.csv",head=T,sep=",")
groupf=unique(file$con.file.name)

file1=file[file$con.file.name==groupf[1],]
file1$HMM.group=paste0("NO.",file1$HMM.group)
filee=merge(file1,anno,by="HMM.group")

###
#library(metafor)
#tdata=data.frame(HMM.group=filee$HMM.group,filee[,4:7])
#tdata1=tdata[tdata$HMM.group=="NO.1",]
#test=meta.MH(case_in_region,case_not,con_in_region,con_not,names=HMM.group,data = tdata1)
###
col_names=c("HMM.group","OR","lower","upper","p.value")
result=data.frame(matrix(NA,15,ncol=5))
names(result)=col_names
group.HMM=unique(filee$HMM.group)
library(meta)
for(i in 1:15){
tdata=data.frame(HMM.group=filee$HMM.group,filee[,4:7])
tdata1=tdata[tdata$HMM.group==group.HMM[i],]
 tdata1$case_not=tdata1$case_in_region+tdata1$case_not
 tdata1$con_not=tdata1$con_in_region+tdata1$con_not
metaor3<-metabin(case_in_region,case_not,con_in_region,con_not,data=tdata1,sm="OR",studlab = HMM.group) 
OR=exp(metaor3$TE.fixed)
upper=exp(metaor3$upper.fixed)
lower=exp(metaor3$lower.fixed)
result[i,]$HMM.group =group.HMM[i]
result[i,]$OR=exp(metaor3$TE.fixed)
result[i,]$upper=exp(metaor3$upper.fixed)
result[i,]$lower=exp(metaor3$lower.fixed)
result[i,]$p.value=metaor3$pval.fixed
}
result=merge(result,anno,by="HMM.group")
result$group="nonsignificant"
 result[result$p.value<0.05,]$group="significant"
 result$OR=log2(result$OR)
 result$lower=log2(result$lower)
 result$upper=log2(result$upper)
 result=result[order(result$num),]
 result1=result[1:9,]
 result2=result[10:15,]
 result=rbind(arrange(result1,group,OR),arrange(result2,group,OR))
 result$num=1:15
ggplot(data=result,aes(y=OR,x=num,group=group,color=group))+geom_point(aes(shape=group),size=3)+
    geom_errorbar(aes(ymin=lower,ymax=upper),width=0.1,color="black")+coord_flip()+
    scale_x_continuous(breaks = result$num,labels = result$description)+theme_light()+ylab("log2(OR)")+geom_hline(yintercept = 0,color="black")+
    scale_color_manual(values = alpha(c('green','red'),0.9))+theme(panel.grid.minor = element_blank())



